4 Dirac-Fock 

4.1 The energy-expression 

Subtracting the rest energy tykP from E, the one-electron Dirac Hamiltonian 

Z 



hi = ca^ ■ p. + mc {Pi - I) 



and the Dirac-Fock Hamiltonian 

i i<j 

The total wave function is the Slater determinant 

* = (iV!)-V2|^i(i)^2(2) . . . Viv(iV)|, 
where the one-electron wave functions 



Then the total energy 



(4.1) 



(4.2) 



(4.3) 



(4.4) 



i<3 



i i<j 



(4.5) 



in terms of the one-electron integrals /, Coulomb integrals J and exchange 
integrals K. 
Eq. (3.129) 



ic 



dr r 



53 



+/ 



f d.{p[c( 



dr 



+Q 



dr r r r 

( dP q q K ^ 
V dr r r r 



P , 1„, P IdP 1 

5=-, 9' = -P'-^ = --, 9 

r r r-^ r ar r 

We see that Inn is independent of m. 
Consider next the general two-electron integral 



(4.6) 
(4.7) 



(■^aV'bl — IV'cV'd) = {ah\v\cd) 
ri2 



dr, / dr2V:(l)V'6(2)— Vc(l)^d(2) 
ri2 



(4.8) 



using the expansion 



1 A; 



E ^ E (-l)^C','(l)C^,(2). 



fe > q=—k 



(4.9) 



Here 



Then the angular parts 



(4.10) 



{Km\Cym') = (-)™+^ \{2j + l)(2j' + 1) 



1/2 



X 



= d {j'm , jm) 



(4.11) 
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and the radial Slater integrals 



/* oo /* oo k 

R\ahcd)= dn dr2^[Pa{l)Pc{l) + Qa{l)Qc{l)] 

Jo Jo 

X [Pb{2)Pd{2) + Qb{2)Qd{2)] 

/•oo 1 

= / -Ykiac, rs) [^6(2)^^(2) + Qb{2)Qd{2)] drs 

Jo ^2 



where 



lfe(ac,r2) =r2^"-^[P„(l)Pc(l) + Qa(l)gc(l)]dri (4.13) 



(4.12) 



{ab\v\cd) = ^{-)'}R''{abcd){a\C^\c){b\C'lg\d). 

kg 

The 3- J symbols fulfil 

\j3m3) = {2js + 1)V2 V {-)-n+h-ms fjl J2 J3 



mi,m2 



(4.14) 



\jimi)\j2m2). 



(4.15) 



In other words, 



C(iij2j3;mim2m3) = (-)-^'i+^'^-"^^(2j3 + l)'/M ] , (4.16) 

\mi m2 —1713 J 



Jl J2 J3 

mi m2 ms 



(_)ji-j2-m3(2_^3 + l)-'/'C(iii2j3;mim2 - ms). 

(4.17) 



The 3-j symbols are tabulated by Rotenberg et al. For their symme- 
tries, see p. 107 of Kvantkemi I. They can be calculated by a simple program. 
In the special case ji = l-, 32 = \i C{l^j;mimsm) is: 



21+1 



21+1 



l-m+^ l l+m+\ 



21+1 



21+1 



mi = m — ms- 



Returning to Eq. (4.11), if the quantum numbers m and m' are fixed, 
only one value q = m — m' survives. 
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Prom 

{ab\v\cd) =^{-yR{ahcd){a\C''q\c){h\C'lq\d), [4.14] 

kq 

{ah\v\ah) = ^{-YR{ahah){a\C';^\a){b\Ctq\h) 

kq 

k 

= ^a''{jama,jbmb)F^{a,b), (4.18) 
k 

where 

a'^UamaJbrrib) = d''{jama,jama)d''{jbmb,jbmb) (4.19) 

and 

F^{a,h)=E!'{ahah) (4.20) 
This was the Coulomb integral. The exchange part 

{ab\v\ha) =^{-yR{abba){a\C^\b){b\C^_q\a) 

kq 

= ^b\jama,jbmb)G\a,b) (4.21) 
k 

with 

b\jm; j'm') = [d''{jm; j'm')f (4.22) 



and 



G^{a,b)=R''{abba). (4.23) 
Using the symmetry rules of the 3-j symbols, 



n 32 .73 \ I 31 33 32 

,mi m2 1713 J XnT-i "rns '7i2> 



J = 31+ 32 + h 



{-Y I I (4.24) 

-mi —m2 —m^j 
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and noting that q is fixed hy q = m — m', 
d''{j'm';jm) = 



-r+^[(2i+i)(2/+i)]^/^ (^'^ n 
-r+i[(2j + i)(2/ + i)]^/M ^' I 
-r+i[(2j + i)(2/+i)]^/^ f^'^ I 

U -1/ 



-yd''{jbmb;jama). 



d'^UamaJbrnb) -- 

For one electron outside a filled shell, j, 



j k j' 
-m q m' 

j' k j 
^m! q —my 

/ k j 
. —m' q mj 

(4.25) 

(4.26) 
(4.27) 



Thus the A; > Coulomb terms vanish in case the potential is spherical. The 
exchange in this case, 



(4.28) 



with 



j k f 



1 -i 

.2 ^ 2/ 



The total energy for a closed-shell atom becomes 

E^ = Y,QaI{A), qA = 2jA + l, 

A 

E'' = Y.{llAiQA-l)F'{A,A) 



(4.29) 



(4.30) 
(4.31) 



A 



1 ^qA-l 



- k>0 ^ V ' 



+ 1 E QAqB{F\A,B)-Y,l^j.kj,G'^{A,B)] (4.32) 



A,B AjtB 
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a. k = part of the intrashell repulsion 

b. 1 for a full shell, for one electron 

c. Exchange of one Ja electron with the filled shell A 

d. k = intershell part (the only Coulomb term) 

e. The exchange of one Ja electron with the average shell B electron 



Example 11 The intrashell interactions of one shell: 

All radial integrals are identical. The angular parts (m m') become 



Case a : A; = 

k = 0^q = 0^m = m'. No self-interaction the exchange term 
lP{jm, jm!) = 0. 

There are six aP{m,,m') in (4.34). Their value equals 1. See LP. Grant, 
Proc. Roy. Soc. Lond. A 262 (1961) 555, Tables 1-3 or use (4.19)-F(4.25). 




(4.33) 




(4.34) 



Note that there are no self-interacting terms, a^{\ |) etc. 
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E.g. 



a^i- -) = cP(- - - -)(f(- - --) 



-16 

-16- 

1 □ 



(2^ + 1) 



3 3 
2 2 



1 1 

'2 2 



Ix 1 



T2 




See KK I, p. 109 



(4.35) 



(Grant, op.cit., Table 2, same result). 6x1 = 6. The first term of E'-' , eq. 
(4.32), gives 



^qiq - 1) = 1 . 4(4 - 1) = 6 ^ 6F°(A A). 



(4.36) 



Case b: A; = 2 



Similarly, 



111111 

+ 7:^ + 



25 25 25 25 25 25 

2 

" ~25 

Symmetries: 



(4.37) 



a'^{jm; j'm') = a^{j'm';jm) = a^{j — m;j' — m) 

= a^{j - m; j'm') = a''{jm; j' - m!) (4.38) 
h^{jm; j'm!) = h^{j'm'; jm) = h^{j - m; f - m!) (4.39) 
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The coefficients are given in Table 4 of Grant. Here, 



2 

25 
8 

25 



+ 



25 



+ 



+ 



+ 



25 



2 

25 



Tfie radial parts R{abab) = R{ahha) = R{aaaa). Adding, 



i:«^-i:'^=(-|)-(|)=-^ 



Ineq. (4.32), 



],(^)(2. + i)ir.,.. -1.4.^.4.1.1 



(4.40) 



(4.41) 



(4.42) 



Grant, Table 5, p. 570 



4.2 The Dirac-Fock equations 

Tlie two-electron energy (4.32) equals 



E 



c 



I 

\ {qa - l)r- 



oo ^ 

r L 



PA{r) + Q\{r) 

oo -j^ 



(4.43) 



-ds 



r> 



+ qsr i — ds 
S^A ^> 



Pl{s) + Q\{s) 
Pl{s) + Ql{s) 



E 

fe>0 



2jA r 



oo 



{2jA + l)Vj^kj^r / ^ Pi{s) + Ql(s) 



r 



ds 



^A^B fc ^ Jo r> 

X [Pa(s)Ps(s) + Qa{s)Qb{s)\ [PA{r)PB{r) + Qyi(r)gB(r 



c. 



d. 



When introducing the variation, AQ, a factor of 4 will appear in each of the 
four terms, from: 

a. A(Q2) = 2QAQ, both AQ(r) and Ag(s) 

b. A(Q2) = 2QAg, = i or 5 = i 
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c. As a. 

d. Q(r) or Q{s), A = i or B = i 



(4.44) 



In terms of the variations AP and AQ, the variation of the two-electron 
energy becomes 



=2qA 



Dr 



PaAPa + QAAQa 



(4.45) 



{(,.-i).f i 



ds 



ds 



Pl{s) + Q\{s) 
PUs) + QUs) 



fe>0 



2jA 2 



{2jA + l)rj^kjAr 



CO „fc 



ri+' 



Pl{s) + Q\{s) 



ds 



poo poo 

J2 ^jAkjB / dr / ds 

[Pa{s)Pb{s) + Qa{s)Qb{s)\ [PB{r)PAPA + QB{r)AQA 



X 



We introduce the auxiliary function 

Yn{A, r)=- {qa - l)Yo{A, A; r) - qBYo{B, B; r) 
QA-ll, 



A:>0 



2jA 2 



{2jA + l)rj^kjMA,A;r), (4.46) 



where 



Yk{ac;r) =r 



oo k 
' < 



4+' 



Pais)Pcis) + Qais)Qcis) 



ds 



[4.13] 



could be interpreted as an effective electronic charge, containing both inner 
and outer screening. 

Orthogonality leads to the Lagrange coefficients e, whose variation 



Ae =A 



- ^ qaeaaNaa - ^ QaQbsabNab 

A A^B 



where 



/•oo 

Nab = / (PaPb + QAQB)dr. 

J o 



(4.47) 



(4.48) 
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Introducing ea = rnc? — Ea, eq. (3.134) gives 



c(P' + -P) + (-mc^ + V-mc^-s)Q = 

r 

c{Q' - ^Q) + (-mc^ -V + mc^ +e)P = 



+ Qx 

- Px 



(4.49) 



e(p2 + Q^) = c\qp' + -QP - PQ' + -PQ 
I r r J 

+y[p2 + q2] 

= c gP' - PQ' + yPg - 2mc2g2 + V[P^ + Q%{A.bO) 



The variation of this one-electron term, 



Asa = QA 



2k 

AQP' + QAP' - APQ' - PAQ' + —{PAQ + QAP) 



-Ac^QAQ + 2F(PAP + QAQ)}dr. 



(4.51) 



Prom AP = AQ = we get by partial integration 



APdr 



1*00 poo 

/ QAP'dr = - / Q' 

Jo Jo 

/ PAQ'dr = - / P'AQdr 
. Jq Jq 



(4.52) 



Hence the variation of the one-electron energy 



Asa = 2qAC | AQaP'a + AP'^Q'a + - {PaAQa + QaAPa) 
-2cQaAQa + ^V{PaAPa + QAAQA)}dr. 



(4.53) 



Adding to this the variation of the two-electron energy, AE'-^ , eq. (4.45), 
and sTibtracting the Lagrange multiplier terms for normalization and orthog- 
onality, we get 
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r IqacI/^QaPa - APaQ'a + —{PaAQa + Qa^Pa) (4.54) 
Jo ^ f 

- 2cQaAQa + WniPAAPA + QaAQa)} 



Yn{A,r) 



-2qA / — [PaAPa + QaAQa 

roo J 

- E E ^^Akjs / -YkiA, B- r)[PBAPA + QbAQ, 

B^A k Jo f 

roo 

-2qASAA / iPAAPA + QAAQA)dr 

Jo 

roo 

-2 V qAQBSAB / {APaPb + AQAQB)dr = 0. 

D^A Jo 



a. 



b. 



d. 



e. 



B^A 

a. one-electron 

b. use (4.46) 

c. (4.45), inter-shell exchange 

d. normalization 

e. orthogonality 

In the last term (e.), a factor of 2 comes from the possibilities "A" = A 
and "5" = A. Dividing (4.54) by 2qAC the coefficient of AQa and -APa 
become 



P'a + —Pa + -[-2c2 + Vn- -Yn{A, r) - £aa]Qa 
r c r 

+ Y 1b[-Y ^jAkjB ^Yk{A, B; r) - -eabK'^a, hb)\ Qb = 
B^A k 



Q'a + —Qa + -[-Vn + -Yn{A, r) + eAA]PA 

ly* ^ 



+ 



1 



1 



Bj^A 



^jAkjB 2^^^^"^' -£ab6{ka, Kb 



Pb = 



(4.55) 

These are the Dirac-Fock equations for a closed-shell atom. We may set 
^Y{A, r) = -Vn + ^Yn{A, r) = ^[Z + Yn{A, r)] (4.56) 
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4.2.1 The relativistic Koopmans theorem 

Multiply (4.55a) by cQa, (4.55b) by —cPa, add, integrate and suppose 
EAB = 0: 



j^" cQa (p'a + ^Pj) - cPa (q'a - "yQA) - 2^Q\ (4.57) 
+ Vn{Pl + Q\) - SAAiPl + Q\) + {Pi + Q\)(--Yn{A,r)) 
+ c 1b{QaQb + PaPb) J2 ^jAkjB -^^k{A, B- r) 



B^A 



dr = 0. 



BAA = £ [cQa [p'a + ^Pa) - cPa (q'a - IQa) - 2c^Q\ + K 

/•OO Q._\ 1 

+ dr{Pl + Q\) +'i^Yo{A, A-r)+Y, -QbYo{B, B; r) 



Ia 

[4.46] 



E \ + l)r,,,,,n(A, A- r) 



k>0 



2jAr 2 



+c QB{QA{r)QB{r) + PA{r)PB{r)) ^JAkjs ^ 



B^A 



I 

Jo 



OO „k 



2cr 



[PA{s)PB{s) + QAis)QB{s)]ds 



[4.13] 



= lA + {qA-l)F\A,A)+YQBF\A,B) 

B^A 

- E + l)r j.kj.F'' {A, A) 



k>0 



2jA 2 



E ^-Y.^j.k,,G\A,B). 



(4.58) 



St^a k 



A comparison with the total energy, (4.32) shows this to be the difference 
E{qA)-E{qA-l) □ 

4.2.2 Multiconfiguration treatment, a simple example 

(C. Briangon, J.P. Desclaux, Phys. Rev. A 15 (1976) 2157) 
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The total energy for a p^{J = 2) configuration. At the LS-hmit, 



1 

1 
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Generally, 



I J = 2) = a|pV) + Vl-a^N^), 
with the total energy 

Et = a^E^^ipV) + (1 - a^)E^y{pp^) 



(4.59) 



(4.60) 



+ - V2aVl - a?R^ {p, P, P, P) . 
4.2.3 " Average-of-configuration" treatment 



(4.61) 



Assume that all \ jm) states of the one-electron states are equally probable 
in the many-electron wave function. 

Example 21 Iodine, bp^ . 

The hole can be in the hp or the 5p shell. They have 4 and 2 mj-states, 
respectively. Therefore the statistical weights are g and |, respectively. 

Example 32 Ground state of uranium, 5/^6d. 
a) /-configuration. 





Z = 3, 


f--j = ^+\- 


7 

" 2' 


2i + 1 = 


8 










5 

" 2' 


2j + 1 = 


6 




config 




Pf 




IP 




E 


weight 


^ = 20 


f • ^ = 120 


6 

1! ■ 


If = 168 


^ = 56 


364 




_ 5 
91 


_ 30 
91 




42 
91 


14 
91 





Note that you have to include the permutation of electrons, n!, in the 
weight calculation. 

b) d-configuration 

1 = 2, d:j = 2 + ^ = ^, 2j + l = 6 

d:j = 2-^ = ^, 2j + l = 4 
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config 


d 


d 


E 


weight 


4 


6 


10 




2 


3 






5 


5 





Altogether eight jj-coupled configurations: 



config 

weight 
weight X 455 


dp 

2 5 
5 ■ 91 

10 


dpf 

2 30 
5 ■ 91 

60 


dfP 

2 42 
5 ■ 91 

84 


^7=^ 

2 14 
5 ■ 91 

28 


config 
weight 
weight X 455 


dp 

3 5 
5 ■ 91 

15 


dpf 

3 30 
5 ■ 91 

90 


dfP 

3 42 
5 ■ 91 

126 


dp 

3 14 
5 ■ 91 

42 



4.3 Numerical solution of the Dirac-Fock equations 

Choosing the variable 

i = Inr 



(4.62) 



and the step length h = At = 0.05, one usually obtains a precision of about 
9 figures (for Et)- About 200 radial grid points, rj, are needed between 
the first point and the practical infinity. All wave functions Pa, Qa and 
potentials lY{A,r) and lY^Bj^AlsY^k^jAkjB^kiA, B; R) are expressed at 
the Vi. 

The integrations outwards are started by the series expansion 



P{r) =r'^{po+pir + p2r^+ ■ ■ ■) 
Q{r) = r'^{qo + qir + q2r^+ ...). 



(4.63) 



After sufficiently many points are obtained, a five-point Adams predictor- 
corrector method or some other numerical method is used. Consider 



y 



dt 



f{y)y{t)+9{t) 

The new, (n + l):th point is obtained by combining the predictor 



(4.64) 



Pn+i =yn + ^ [I901y; - 2774y;_i + 2616y;_2 - 1274y;_3 + 251y;_4] 

(4.65) 
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and the corrector 



Cn + l=yn + ^ [251^+1 + 646y; - 264y'^_^ + W6y'„_^ - 19y'„_3] 

(4.66) 

to 

Vn+i = ^ [473c„+i + 29pn+i] + 0{h^) (4.67) 

The outwards integration becomes unstable after the (effective, radial) 
potential, V{r) = —^Y{A,r) crosses the eigenvalue, £ 




Figure 4.1: The instability of the integration depicted 



Similarly, the inwards integration becomes unstable after this point. The 
solution is to fit the two integrations at this point. 

In "Method 2" one solves both the original inhomogeneous differential 
equations and the homogeneous ones with YIb^a ~ ^- Then every linear 
combination of the two solutions is a solution. 

The coefficients of the homogeneous solution, a, follow from 

Pi(ri) + ai^Ptin) = P^(ri) + a^Pkin) 
QL(ri) + aLQ^(ri) = Qk(n) + aRQ^(ri), 

at the matching point, rj. Here i=inhomogeneous, h=homogeneous, L="left" 
(r < n), R=" right" (r > n). 

This solution fulfills the boundary conditions but may not be normal- 
ized. Normalization by a constant is not allowed because the equation is 
inhomogeneous =^ the eigenvalue of eaa is chosen to obtain normalization. 
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Actually it suffices if we converge towards a normalized solution. Let 
£ = eaa- Then 

+ Q^)dr + 25e + Q^)dr = 1. (4.69) 

By first differentiating the original Dirac-Fock equations, 



dAdeJ^'^de^^^'^' de- 2c^ 



(4.70) 



these derivatives can be solved by the same subroutines as P and Q them- 
selves. Then (4.69) 



2r(^'^S+0Aff)d'-' 



In the old Dirac-Fock One-Center Expansion (DF-OCE) method, four 
numerical methods could be used: 



1. • Calculate only ip\, put Pl = -Pr + "r^r- 

• Fit Q by changing e. 

• Finally normalize by force (inner shells) or by using dP/de (molecules). 

2. Previous method. Normalize by adding to P some amount of 5e x 
(dP/de). If the number of nodes requires it, change e. The continuity 
is guaranteed. 

• Fit e by (4.71) 

3. Use both ol and or. Use (4.71). For wrong number of nodes, 

4. change e to get the right number of nodes. 



4.3.1 Specific features of the DF-OCE method 

• Each MO consists of a dominant AO and other AOs. Example: For 
the Fy MO of CH4, 2p3/2 + 4d + 4/. 

• Use for the dominant component Method 1 and for the others Method 
3 or, if necessary. Method 4. 

• Then the dominant component fixes £ (which is not changed in Method 

3 for the other components). 
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After treating the entire MO, set the common 

^MO = Yl ^(AO)£AO 
AO 

(SOLDMO, hnes 175-177) 
Then correct it by setting 



£mo — £mo + 



Pi 



(SOLDMO, hnes 178, 180) 



(4.72) 



(4.73) 



• After treating the entire MO, normahze it by the coefficient Ni] 
(SOLDMO, hnes 157, 204-211) 



Some results: 
Desclaux and Pyykko (1974) 

D&P (1976) 

D&P (1976) 

D&P (1977a) 

D&P (1977b) 

D&P (1978) 

P (1979b) 

P (1979c) 
[See Table 7.3 of RTAM I.] 



the bond-length contraction 
the Ag/Au difference 
p bonding in TIH. 
TiH4-(104)H4. Trends. 
(114)H4 > PbH4. 

MoHe/WHe- Re,De,k2; oxidation state. 
MH+, MH2, Groups 2 and 12. 
LnH, AnH. 



The iterations can be speeded up (by roughly a factor of 2) by using a 
convergence factor a. 

Let Pn,Qn and Pn,Qn be the initial and final values during the n:th 
iteration. Choose 



f PL+i{r) = P[aPl + {l-a)P^] 
\Ql+iir) = P[aQi + {l-a)Ql] 
Here /? is a normalization factor and a is restricted to 

0.1 <a< 0.9 

and is calculated every time from the change 

ARn = max[P^(r) - P^(r), Ql{r) - Qi{r)]. 



(4.74) 
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